% Estimation of the logp1 with Area dependce models
% s3: 
% REML estimator is presented
% 3/1 version

%% Housekeeping
clc; clear all; close all;
workpath00 = pwd;

addpath('toolbox_xt');
addpath('toolbox_plot');

%% load data
load('data_all.mat');

%% Two panel variables
% xtflag and xtval: combofips
% xtflag2 and xtval2: year
% xtflag3 and xtval3: combofips/year
xtflag3 = xtflag + xtflag2/10000;
xtval3 = unique(xtflag3);

nobs = size(YY,1);
ncity = length(xtval);
nyear = 6;

%% Data for the esitmation
% sel_est = true(nobs,1);  %include all

[~,~,info] = xt_reg(YY, ones(nobs,1), xtflag3, xtval3);
% sel_est = (info.n>=5);
sel_est = true(size(YY,1),1);

est.YY  = YY(sel_est,:);
est.xtflag = xtflag(sel_est,:);
est.xtval  = unique(est.xtflag);
est.xtflag2 = xtflag2(sel_est,:);
est.xtval2  = unique(est.xtflag2);
est.xtflag3 = xtflag3(sel_est,:);
est.xtval3  = unique(est.xtflag3);
est.nobs    = size(xtflag,1);

% Other controls
YD = [year_dummy2, year_dummy3, year_dummy4, year_dummy5, year_dummy6];

% estimation
est.YD = YD(sel_est,:);
est.XX = XX(sel_est,:);
est.logDC = logDC(sel_est,:);
est.DD = DD(sel_est,:);

est.ZZ = ZZ(sel_est,:);
[~,temp_ZZ2] = xt_reg(est.ZZ(:,2), ones(size(est.ZZ,1),1), est.xtflag, est.xtval);
est.ZZ2 = [ones(size(temp_ZZ2,1),1), temp_ZZ2];

estinfo.islog = 2;
estinfo.FitMethod = 'REML';
% % unshrunken
% est_us_m2 = estimator_us_p1(est, estinfo);
% est_us_m2.model.MSE
% 
% est_us_m5 = estimator_us_p4(est, estinfo); %almost impossible to do this esitmation...
% est_us_m5.model.MSE
% 
% % shrinkage
% est_s_m2 = estimator_s_p1(est, estinfo);
% est_s_m5 = estimator_s_p4(est, estinfo); %almost impossible to do this esitmation...

% Estimation of the models
est_s1 = [];
est_s2 = estimator_s_A_p1s3(est, estinfo); %
% est_s3 = estimator_s_A_p2(est, estinfo); %
est_s3 = [];
% est_s4 = estimator_s_A_p3(est, estinfo); %
est_s4 = [];
% est_s5 = estimator_s_A_p4(est, estinfo); %
est_s5 = [];
est_s6 = estimator_s_A_p4s3(est, estinfo); %
% est_s7 = estimator_s_A_p4sr(est, estinfo); % random in size

estinfo.FitMethod = 'ML';
est_s7 = estimator_s_A_p4s3(est, estinfo); %
% est_s7 = '';

% save
save('results_full_logp1_A_s3_REML.mat', 'est_s1', 'est_s2', 'est_s3', 'est_s4', 'est_s5','est_s6', 'est_s7');



